% genhardcore(n, lmd, K)
%
% simuliert das Allgemeine Hardcore Modell mittels MCMC
%
% IN:
%       n   - Laenge bzw. Breite einer Konfiguration
%       lmd - Lambda, steuert die Anzahl schwarzer Pixel
%       K   - Anzahl der Iterationsschritte
%
% Beispielaufruf:
%
%       >> genhardcore(8,1,10000)

function [] = genhardcore(n, lmd, K)
    X = zeros(n,n);             % G?ltige Anfangskonfiguration
    N = sum(sum(X));            % Anzahl schwarzer Pixel

    for k = 1:K
        i = ceil(rand(1,2)*n);  % Zuf?lliger Pixel

        if (isPixelFree(X,i))
            if X(i(1),i(2)) == 0
                if rand(1) <= lmd
                    X(i(1),i(2)) = 1;
                end
            else
                if rand(1) <= 1/lmd
                    X(i(1),i(2)) = 0;
                end
            end
        end
                        
        N = [N sum(sum(X))];    % Anzahl schwarzer Pixel
    end

    
    % Plotte Ergebnisse
    
    figure()
    
    subplot(1,4,1)
    colormap(hot(256))
    image(255*(1-X));
    title('X_K')
    
    subplot(1,4,2)
    plot(N)
    xlim([0 K]);
    ylim([0 n^2/2]);
    title('Anzahl schwarzer Pixel')
    
    subplot(1,4,3)
    plot(cumsum(N)./[1:length(N)]);
    xlim([0 K]);
    ylim([0 n^2/2]);
    title('gleitendes Mittel der Anzahl schwarzer Pixel')
    
    subplot(1,4,4)
    hist(N(500:10000),[min(N):max(N)])
    xlim([0 n^2/2]);
    title('Histogramm')
    
    mean(N)
    mean(N(1000:end)) % mit Burn-In Phase
    
end

% Liegen in der Nachbarschaft von Pixel p in der Konfiguration X schwarze
% Pixel ?
function free = isPixelFree(X, p)
    n = size(X,1);
    m = 0;
    if (p(1)>1)
        m = m + X(p(1)-1,p(2));
    end
    if (p(1)<n)
        m = m + X(p(1)+1,p(2));
    end
    if (p(2)>1)
        m = m + X(p(1),p(2)-1);
    end
    if (p(2)<n)
        m = m + X(p(1),p(2)+1);
    end
    free = (m == 0);
end